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ABSTRACT 


Active  sensing  in  structural  health  monitoring  (SHM)  refers  to  injecting  (user-defined)  energy  into 
the  system  in  order  to  actively  probe  its  response  to  the  induced  dynamics  as  a  means  of  detecting 
whether  damage  may  be  present  in  the  system.  A  number  of  researchers  have  shown  that  active 
sensing  with  guided  ultrasonic  waves  (GUWs)  can  be  a  powerful  approach  to  take,  as  GUWs,  when 
launched  and  detected  in  conjunction  with  macro-fiber  composite  (MFC)  patches,  can  retain  the 
wide  area  coverage  capability  of  lower  frequency  (vibration)  methods  while  significantly  enhancing 
sensitivity  to  flaws  (cracks,  corrosion,  etc.)  because  of  the  small  interrogating  wavelengths  used. 

This  project,  led  by  PI  Michael  Todd  and  Co-PI  Francesco  Lanza  di  Scalea,  considered  several  new 
concepts  in  guided  GUWs:  (i)  optimized  passive  GUWs,  (ii)  quantitative  active  GUWs,  and  (iii)  a 
generalized  insonification  (diffuse  field)  approach  rooted  in  data-based  modeling  and  pattern 
recognition.  These  concepts  arc  tested  on  metallic  test  articles  to  detect  a  variety  of  defects, 
including  impact  damage,  simulated  cracks  (notching),  corrosion,  and  bolted  joint  preload  loss. 

1.  INTRODUCTION 

The  use  of  original  design  solutions  and  innovative  materials  in  response  to  the  evolving 
demands  of  the  industry  introduces  challenges.  In  the  naval  industry,  for  example,  there  has  been 
the  development  of  high-speed  and  high-performance  vessels.  The  aluminum  material  being 
considered  for  the  newest  Navy  platforms  is  susceptible  to  stress  corrosion  cracking,  fatigue 
cracking,  or  sudden  failure  due  to  the  impact  of  waves  or  small  icebergs,  called  bergy  bits. 
Bolted  connections  in  a  number  of  structural  components  are  also  of  interest. 

Structural  health  monitoring  solutions  able  to  detect  and  quantify  structural  defects  could  be 
beneficial  in  these  innovative  naval  structures  providing  information  on  the  condition  of  the 
structure  in  real-time  (Condition-Based  Monitoring).  Overall,  the  approach  is  an  "active  and 
passive  sensing"  approach  whereby  elastic  waves  are  created  and  propagated  through  the  hull, 
and  features  derived  from  the  waves'  interaction  with  the  hull  arc  mined  for  evidence  of  hull 
condition  changes.  While  the  passive  approach  is  a  form  of  acoustic  emission  whereby  growing 
cracks  are  detected  in  real  time,  the  active  sensing  approach  is  meant  to  detect  and  quantify  both 
growing  and  stable  cracks.  Furthermore,  a  third  form  of  elastic  wave  probing  (in  the  active 
mode)  that  takes  advantage  of  the  diffuse,  reverberant  field  (rather  than  explicit  mode  probes)  is 
explored  as  a  means  of  detecting  and  classifying  defects. 

In  the  passive  approach,  growing  damage  or  sudden  impacts  can  be  located  by  monitoring 
acoustic  emissions  in  a  real-time  mode  (Mai  et  al  2005,  Banerjee  et  al  2005,  Prosser  et  al  1995). 
Traditionally,  damage  or  impact  location  is  based  on  time-of-flight  triangulation  of  wave 
measurements  taken  at  multiple  receiving  points.  One  problem  with  time-of-flight  triangulation 
is  the  dispersive  nature  of  the  flexural  mode  which  makes  it  difficult  to  identify  an  accurate 
arrival  time  from,  for  example,  the  conventional  first-threshold  crossing.  This  concern  remains 
despite  the  proposed  improvements  over  the  first  threshold  crossing,  such  as  cross-correlation 
techniques  for  extracting  the  phase  of  a  single  frequency  component  (Ziola  and  Gorman  1991). 
The  measurement  of  the  less  dispersive  extensional  mode  can  theoretically  help  identifying  a 
correct  arrival  time  (Gorman  and  Ziola  1991).  However,  the  extensional  mode  can  be  highly 
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attenuated  in  composites  or  generally  large  structures,  in  which  cases  triangulation  of  the  flexural 
mode  often  remains  the  only  realistic  option.  Another  reason  for  preferential  use  of  the 
(fundamental)  flexural  mode  is  the  lower  phase  velocity  compared  to  the  (fundamental) 
cxtensional  mode  at  the  same  frequency.  This  results  in  smaller  wavelengths  thus  higher 
sensitivity  to  small  damage.  The  second,  more  severe  problem  with  time-of-flight  triangulation  is 
the  required  knowledge  of  the  wave  velocity  in  the  test  material.  The  velocity  is  required  to 
translate  arrival  time  measurements  into  source  location.  This  requirement  is  a  fundamental 
limitation  when  monitoring  anisotropic  or  geometrically-complex  structures  (c.g.  aerospace 
panels),  where  the  velocity  is  dependent  upon  propagation  direction.  Model-based  approaches 
(Chang  and  Markmiller  2006,  Park  and  Chang  2005,  Seydcl  and  Chang  2001)  have  thus  been 
required  for  impact  location  and  identification  in  simple  composite  structures.  In  these  methods 
the  location  and  type  of  impact  in  a  model  are  iteratively  changed  until  the  predicted  responses 
match  the  measured  responses.  For  damage/impact  location  in  more  complex  structures,  where 
accurate  model-based  predictions  may  not  be  possible,  artificial  neural  networks  have  been 
proposed  (Staszewski  et  al  2000,  Jones  et  al  1997,  Schindler  et  al  1995).  However,  neural 
networks  require  an  extensive  number  of  training  observations  prior  to  deployment,  and  they 
cannot  provide  deterministic  locations.  However  all  the  works  previously  mentioned  has  focused 
on  impact  tests  performed  in  a  low  velocity  regime.  In  fact  the  actual  impacts  were  generated 
using  hand-held  hammer,  falling  weight  or  just  pencil  lead  breaks  (Hsu-Niclsen).  What  arc 
missing  are  investigations  of  impacts  in  a  more  realistic  velocity  regime.  This  paper  presents 
experimental  investigation  of  high  velocity  ice  impacts  location  based  on  piezoelectric  rosettes. 
Ice  projectiles,  were  projected  at  velocity  around  30  m/s  using  a  gas-gun  onto  a  6061-T6 
aluminum  plate  and  a  [0/+45/90/-45]s  woven  composite  plate.  The  rosettes  exploit  the  highly 
directional  response  of  Macro  Fiber  Composite  (MFC)  rectangular  patches  to  determine  the 
direction  of  an  incoming  wave  independently  of  the  wave  velocity  in  the  medium.  MFC  patches, 
originally  developed  at  NASA  Langley  Research  Center  for  low-frequency  structural  control 
(Wilkie  et  al.  2000,  Sodano  et  al.  2004),  are  recently  being  used  for  guided-wave  transduction  in 
aircraft  structures  owing  to  their  superior  flexibility  and  durability  compared  to  monolithic  PZTs 
(Lanza  di  Scalca  et  al.  2007a,  Matt  and  Lanza  di  Scalca  2007,  Salamonc  et.).  Since  no  wave 
velocity  information  or  complex  modeling  is  required,  the  proposed  technique  can  improve 
source  location  in  anisotropic  and  gcometrically-complex  structures.Thc  general  idea  of  using 
rosettes  for  the  detection  of  ultrasonic  waves  is  not  new.  However,  the  use  of  piezoelectric 
rosettes  for  extracting  the  direction  of  incoming  waves,  thus  allowing  real-time  source  location, 
is  novel.  Recent  works  employing  rosettes  of  Bragg  grating  fiber-optic  strain  sensors  (Thursby  et 
al  2004  and  2006)  for  wave  direction  determination.  However,  current  technology  for  Bragg 
grating  detection  of  ultrasound  requires  manual  tuning  of  the  laser  to  a  high-sensitivity 
wavelength  specific  to  each  grating.  Consequently,  ultrasonic  measurements  cannot  be  taken 
from  all  three  gratings  of  the  rosette  simultaneously,  which  limits  the  technique  to  locating 
existing  defects  only,  rather  than  growing  damage  or  impacts  in  real  time.  In  addition,  the 
periodicity  of  the  Bragg  grating  response  requires  sophisticated  data  reduction  steps,  such  as 
genetic  algorithms  (Thursby  et  al  2004),  to  identify  univocally  the  wave  incoming  angle.  These 
issues  are  avoided  with  piezoelectric,  rather  than  optic,  wave  sensing.  In  a  previous  study  (Matt 
and  Lanza)  the  performance  of  the  rosettes  for  source  location  was  validated  through  pcncil-lcad 
breaks  performed  on  an  aluminum  plate,  an  anisotropic  CFRP  laminate,  and  a  complex  CFRP- 
honeycomb  sandwich  panel.  In  this  work,  experimental  results  will  be  shown  to  demonstrate  the 
applicability  of  this  approach  to  a  more  realist  impacts  as  ice  projectiles. 
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In  the  active  approach,  guided  waves  are  injected  in  a  defected  region,  thereby  inducing  complex 
scattering.  By  studying  the  characteristics  of  the  reflected  and  transmitted  waves,  quantitative 
information  on  the  defects  is  obtained.  Quantitative  defect  information  can  be  extracted  by  using  a 
hybrid  formulation  wherein  the  finite  clement  method  is  used  to  model  small  regions  near  the  defect 
whereas  regions  away  from  the  defect  arc  modeled  using  a  suitable  set  of  wave  functions.  Dong  et 
al  (2004)  developed  a  global-local  finite  element  formulation  for  modeling  axisymmetrie  scattering 
of  a  steady,  compressive,  incident  elastic  wave  in  a  homogeneous,  isotropic  host  medium  with  an 
axisymmetrie  inclusion.  The  hybrid  method  was  also  used  (Al-Nassar  et  al  1991, 
Rattanawangcharoen  et  al  1997)  to  model  defects  in  a  plate  and  a  cylinder,  respectively.  The 
method  has  been  applied  to  model  wave  interaction  with  defected  lap-shear  joints  (Chang  and  Mai 
1995),  as  well  as  notches  and  rivet-hole  cracks  in  plates  (Mai  and  Chang  2000,  1999).  Doyle 
(Gopalakrishnan  and  Doyle  1 995,  Doyle  1 995)  employed  a  similar  global-local  scheme  wherein  the 
global  field  was  modeled  by  using  spectral  elements. 

Many  components  arc  complex  in  either  their  geometry  (e.g.  multilayers,  tapered  thickness,  etc.)  or 
their  material  properties  (e.g.  anisotropic).  In  this  case  theoretical  wave  solutions  for  the  global 
portion  arc  either  nonexistent  or  hard  to  determine.  The  Semi  Analytical  Finite  Element  (SAFE) 
method  can  help  handle  these  cases  because  of  its  ability  to  extract  modal  solutions  of  complex 
structures  in  a  computationally  efficient  manner  (Bartoli  et  al  2006,  Hayashi  et  al  2003,  Treyssedc 
2008,  Castaings  and  Lowe  2008,  Hayashi  et  al  2006,  Al-Qahtani  et  al  2005,  Lanza  di  Sealca  et  al 
2007,  Matt  et  al  2005).  Sabra  (Sabra  et  al  2008)  demonstrated  the  application  of  the  hybrid 
formulation  to  the  detection  of  holes  in  aluminum  plates.  The  present  paper  extends  the  global-local 
approach  to  model  notches  in  aluminum  plates  and  delamination-like  defects  in  composite  panels. 

This  previous  active  approach  inherently  relies  upon  eharacteristies  of  specific  GUW  modes  and 
their  properties,  and  the  GUWs’  interaction  with  the  structure.  When  coupled  with  SAFE 
models,  very  specific  features  may  be  identified  that  can  indicate  specific  flaws.  However, 
structural  geometries  such  as  highly  curved  regions,  bolted  interfaces,  etc.,  pose  significant 
modeling  and  testing  challenges  for  GUWs,  resulting  in  difficulty  in  classifying  size  or  type  of 
damage  due  to  the  inherent  simplicity  of  the  actively  imparted  excitation  signal,  which  is  usually 
a  modulated,  narrow-band,  short-wave  pulse.  Instead,  some  researchers  have  attempted  to 
employ  bulk  msonification,  where  an  ultrasonic  source  is  excited,  and  the  resultant  long-time,  or 
diffuse,  wave  field  is  examined  to  identify  structural  changes  (Michaels  and  Michaels,  2005). 
This  method  is  investigated  as  well  to  augment  the  standard  guided  wave  method  for  structures 
with  complex  boundary  conditions  or  geometries  that  make  tracking  and  analysis  of  a  single 
propagating  mode  difficult  or  impossible.  This  approach  relies  on  data-driven  (typically,  time 
scries)  models  and  change  detection  approaches  rooted  in  statistical  pattern  recognition. 

Given  these  varied  approaches,  one  might  envision  a  possible  general  taxonomy  for  damage 
detection  strategies  in  the  ultrasonic  domain  as  shown  in  Figure  1.  SAFE  (physics-based) 
modeling  is  critical  for  the  middle  column,  and  data-based  modeling  is  critical  for  the  right-hand 
column.  This  project  presents  studies  that  consider  new  ways  how  the  structural  health 
monitoring  (SHM)  “online”  approaches  can  work  for  metallic  structure  defect  detection, 
localization,  and  classification  and  form  a  complementary  capability  for  Navy  applications. 
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Figure  1  Ultrasonic  domain  diagnostic  taxonomy. 


2.  PASSIVE  APPROACH 


2.1.  Response  Of  Micro  Fiber  Composite  To  Flexural  Lamb  Waves 


Figure  2  MFC  geometry. 


Consider  a  thin  rectangular  MFC  transducer 
(type  P2)  acting  as  a  sensor  and  having  an 
effective  length,  /,  width,  b,  and  thickness,  /, 
along  the  coordinates  1,  2  and  3,  respectively 
(Figure  2).  The  Cartesian  reference  system  ( x , 
y,  z)  is  aligned  with  the  transducer’s 
geometrical  axes  (1,2,  3).  This  transducer  is 
bonded  to  the  upper  surface  of  an  isotropic 
plate  of  thickness  2d  and  subjected  to  a 
harmonic  strain  field  associated  to 
antisymmetric  (flexural)  Lamb  waves 
propagating  in  the  plane  (x  ’,  z)  along  direction 
x  This  wave  mode  is  commonly  used  for  the 
detection  and  location  of  damage  or  impacts. 


i 

2d 

T~ 
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as  diseussed  in  the  introduction.  The  wave  propagation  direction  x  ’  forms  an  angle  0  with  the 
lengthwise  direction  of  the  sensor.  The  origin  of  the  thiekness  coordinate,  z  -  0,  is  at  the  mid¬ 
plane  of  the  plate. 

To  derive  the  MFC  sensor  voltage  response,  in  a  previous  work  Matt  et  al.  modeled  it  as  an 
orthotropic  lamina  in  plane  stress  conditions,  and,  considering  the  eleetrieal  boundary  conditions 
of  the  MFC  sensor  to  be  an  open  cireuit,  the  following  expression  was  obtained,  ignoring  shear 
lag  effects: 


-<(«*+-) 

V  =  iSe,,e  2  (1) 

where  £x.x ,  is  the  amplitude  spectrum  of  the  surface  strain  defined  as: 


”*2 


tanh  rd  — 


2  rs 

kr+72 


-  ■  tanh  sd 


and  S  is  the  frequency-dependent  sensitivity  faetor  that,  for  antisymmetric  waves  is: 


4// 
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■sin 


kb  sin  9 


sin  I 
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(  kl cos  9  ^ 
> 
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(3) 


H  is  a  frequency-independent  constant  depending  on  the  in-plane  orthotropic  elastic  constants  ( 
E, ,  vtj ),  the  elcctro-mcchanical  properties  ( ean  ,  dy ,  dn ),  and  thickness  of  the  sensor  ( t ),  as  well  as 

on  the  wave  propagation  direction  relative  to  the  sensor  geometrical  axes: 


//  = 


+rf„v1!£!)  cos2  6  4-  (d3Xvl2E2  +  dnE2  )sin^  0  J 
[O  —  ^2i^i2  )£?33  —  (5^3 1  +  2<J^dnvnE2  +  dnE2 )] 


(4) 


The  terms  dn  d ^  and  are  the  piezoelectric  constants  and  dielectric  permittivity,  respectively. 
In  this  expression,  obtained  under  the  assumption  of  plain  waves  (ev,y.  =  0),  the  parameters  r  and 
s  are  defined  as: 


where  (O  is  the  frequency,  k  is  the  wavenumber,  and  ci  and  cj  are  the  bulk  longitudinal  and  shear 
velocities  in  the  plate,  respectively.  Matt.  Et  al.  demonstrated  that  the  sensor  behaves  as  an  ideal 
strain  deteetor  when  it  is  small  compared  to  the  wavelength  (i.e.  at  lower  frequencies).  At  higher 
frequencies,  where  the  wavelength  is  comparable  to  the  sensor’s  dimensions,  tuning  effects, 
govern  the  response.  It  can  be  demonstrated  that  forA>/  (sensor  smaller  than  wavelength),  the 

sensitivity  can  be  decomposed  into  a  longitudinal  sensitivity  faetor,  Sx  ,  and  a  transverse 
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sensitivity  factor,  S2 ,  as: 

S  ~  S{  cos2  d  +  S2  sin'  0  (6) 


where  6  is  the  usual  wave  propagation  direction  relative  to  the  sensor’s  lengthwise  direction.  The 
longitudinal  and  transverse  sensitivity  factors  are  defined  for  waves  propagating  along  the 
sensor’s  lengthwise  and  widthwise  directions,  respectively.  Thus: 


St=S 


S2=S 


2'(«+42v12£>n(%) 


Ik  [O' -  v2,v,  2  K  -  «/•,  +  2d,tdnVnE2  +  <£2 )] 
’HdJnE1+dnEl)  sm(^2) 

^  '  M[(l  - v21vl2  K  -  («/*£,  +  2d»dnvnE2  +  d:2E2 )] 


(7) 


By  substituting  Eq.  (6)  into  Eq.  (1)  and  transforming  the  amplitude  spectrum  of  the  surface  strain 
£x.x,  into  the  lengthwise  and  widthwise  axes  of  the  sensor  using  the  following  strain 
transformation: 


cos2  6 
£n=£x.x.sm20 


(8) 


the  sensor’s  response  becomes: 

V  ~  S\  £w  +  S2£22  (9) 


where  the  wave  strain  components  along  the  longitudinal  and  the  transverse  sensor  directions  are 
explicitly  indicated.  Eq.  (9)  makes  the  MFC  response  formally  equivalent  to  that  of  conventional 
electrical  resistance  strain  gages.  This  equation  also  emphasizes  that  the  response  is  highly 
dependent  on  the  wave  propagation  direction,  and  thus  significant  directivity  behavior  is 
expected  at  certain  frequencies. 

2.2.  MFC  Rosettes  For  Wave  Source  Location 

Let  us  consider  three  rectangular  MFC  sensors  A ,  B  and  C  arranged  in  an  arbitrary  rosette 
configuration  and  subjected  to  antisymmetric  waves  as  shown  in  Fig.  2a.  Considering  waves 
such  that  A  >  /  and  using  Eq.  (9),  the  sensors  responses  can  be  written  as: 

VA  =Sxe*+S2e*2 

V  —  ■S’jfn  +  S2£  22  (10) 

V  =  [+ S2£22 
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where and  £22  O'  =A,  B,  C)  indicate,  respectively,  the  strain  components  in  the  longitudinal  and 

in  the  transverse  direction  of  the  i-th  sensor,  and  all  sensors  are  assumed  equal  in  eleetro- 
meehanieal  and  geometrical  properties. 


Figure  3  MFC  rosette  geometry. 


Dividing  Eq.  (10)  by  the  longitudinal  sensitivity  factor  yields 
y‘ 

—  =  £[x+Kt£ 22  for  i  =  A,  B,  C 
where  Kt  is  a  transverse  sensitivity  ratio  defined  as: 

S2  ^  (^31^12^2  +  ^32^2  )S*n  (^2) 

S\  ^(^31£l+t/32Vl2^2)sin(%) 

Let  us  assume  that  the  lengthwise  direction  of  the  z-th  sensor  is  oriented  at  an  angle  a,  from  the 
global  coordinate  axis  *.  In  analogy  with  eleetrieal  resistance  strain  gages,  given  the 
measurement  of  three  independent  strain  from  the  response  of  the  three  sensors  it  is  possible  to 
calculate  the  Cartesian  strain  components  and  their  principal  strain  angle  (p.  The  principal  strain 
angle  of  the  wave,  <p ,  measured  from  the  global  axis  x  (Figure  3a),  can  be  computed  as: 


(ID 


(12) 


tan  2  <p  =  —  (13) 

P  —  P 

xx  yy 

Evaluation  of  the  source  location  in  a  plane  is  readily  achieved  by  the  intersection  of  the 
principal  directions  determined  by  two  rosettes  (Figure  3b).  The  principal  angles  calculated  from 
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the  two  rosettes,  0,  and  02 ,  define  the  equations  of  the  straight  lines  connecting  the  wave  source 
to  each  of  the  rosette  centroids  according  to: 

y SOURCE  ~  (x SOURCE  ~  X\  )  fa\  Tl  (14) 

y SOURCE  ~  (X SOURCE  ~  *2  )  fa  T’ 

where  ( x source  *  y source  )  represent  the  coordinates  of  the  wave  source  in  the  (.v,  y)  Cartesian 
system,  and  (.y,,j>,),  (x2,y2)  represent  the  coordinates  of  the  centroids  of  the  two  rosettes, 

respectively.  The  linear  system  of  equations  can  finally  be  solved  for  the  coordinates  of  the  wave 

source: 


_  >'2  ~  +  xi tan  fa  ~  x2 tan  fa 

X SOURCE  -  "  7  ■  7 

tan  0,  -  tan  0, 

y SOURCE  =  (X SOURCE  ~  Xl  )  *an  fa  Tl 

2.3.  Experiments 


(15) 


In  a  previous  work,  Lanza  di  Scalca  et  al.,  in  collaboration  with  Dr.  Hyonny  Kim,  studied  impact 
location  by  using  two  rosettes  attached  to  isotropic  and  CFRP  anisotropic  specimens,  in  a  low 
velocity  regime  (pencil-lead  breaks).  To  demonstrate  the  potential  of  this  technique,  source 
location  was  attempted  in  the  high  velocity  regime.  A  nitrogen  gas  cannon  (Figure  4),  was  used 
to  shoot  simulated  hail  iee  (SHI)  spheres.  The  monolithic  iee  spheres  had  a  diameter  of  38  mm 
and  were  cast  in  a  split  spherical  mold  over  1  filling  session. 


Figure  4  Gas  gun  for 
high-velocitv  impacts. 


The  specimen  used  for  the  validation  test,  shown  in,  consist  of  the 
0.91  m  x  0.91  m  *  2.5  mm  6061-T6  aluminum  plate.  Two,  three- 
element  P-2  MFC  rosettes,  each  arranged  in  a  delta  configuration 
(120°  between  adjacent  sensors  -  Figure  6),  were  surface-bonded  to 
the  specimen.  During  testing,  the  Fast-Fourier  Transform  was 
performed  for  the  time  history  of  each  sensor.  The  amplitude  was 
computed  in  a  narrow  frequency  band,  near  the  pointA  =  /.  For 
anisotropic  structures,  the  frequency  at  which  A  =  /  is  clearly 
direction  dependent.  Therefore  the  response  magnitude  was 
computed  at  the  maximum  frequency  such  that  optimal  directivity 
response  was  ensured  while  meeting  the  requirement  A  >  /  for  all 
sensors  within  a  rosette. 

A  typical  signal  captured  from  impact  tests  on  the  6061-T6 
aluminum  plate  is  shown  in  Figure  5.  The  flexural  anti  symmetric 
Ao  mode  component  of  this  signal  is  clearly  identified.  This  is 
expected  since  the  impact  produces  a  large  out-of-plane  flexural 
motion. 
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Figure  5  Aluminum  plate  used  in  impact  r 

°  r  r  Figure  6  Typical  recorded  trace  from  impact. 


The  estimated  impact  speed  was  =  30  m/s.  The  estimated  impact  locations  results  are  shown  in 
Figure  7  (right).  In  particular,  Figure  7(left)  shows  a  single  frame  extracted  from  a  high-speed 
film  taken  during  the  impact  test.  As  seen  in  the  picture,  the  iee  sphere  locally  erashes  on  the 
plate  center.  This  aetual  impaet  location  is  marked  in  Figure  7(right)  as  a  red  circle.  The 
predicted  impaet  location,  using  the  proposed  method,  is  represented  on  the  figure  by  a  green 
eirele.  An  RMS  error  function  was  defined  to  provide  a  measure  of  the  rosettes  performance;  the 
RMS  error  is  expressed  as: 


RMS  =  - 
n 


(16) 


where  n  is  the  number  of  the  test  impacts,  {xai,ya,i)  are  the  coordinates  of  the  actual  impaets  on 
the  aluminum  plate  in  millimeters  and  (xpi,ypi)  are  the  coordinates  of  the  predicted  impact 
locations  in  millimeters.  This  form  of  the  RMS  error  function  gives  the  radial  error  (i.e.,  the 
vector  length  from  the  actual  to  the  computed  impact  location).  The  estimated  RMS  error  for  this 
test  was  1 8  radial  mm. 


Length[mm] 

Figure  7  (left)  High-speed  imagery  during  impact  test,  and  (right)  comparison  of  predicted 
and  actual  impact  location. 
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3.  ACTIVE  APPROACH:  EXPERIMENTS  AND  MODELING  USING 
GUIDED  MODES 


The  geometry  of  the  problem  at  hand  is 
depicted  in  Figure  8,  where  the  material  inside 
the  bounded  surface  has  properties  different 
from  the  outside  due  to  the  presence  of  local 
irregularities  as  a  defect.  In  order  to  best 
illustrate  this  hybrid  method,  assume  that  a 
plane  time-harmonic  wave  is  incident  on  the 
FEM  boundary  in  Figure  8.  The  region  inside 
the  boundary  is  meshed  and  analyzed  using 
standard  finite  clement  analysis;  the  region  outside  the  boundary  is  represented  by  the  linear 
combination  of  a  set  of  global  wave  functions,  here  calculated  from  the  SAFE  method.  The 
continuity  of  displacements  and  tractions  on  the  mesh  boundary  must  be  satisfied  through  proper 
choice  of  the  amplitudes  of  the  wave  functions. 

3.1.  Semi-Analytical  Finite  Element  Method 


c 

Incident  wave 

)efect 

•  Ml 

7? 

Transmitted  waves 

Reflected  waves 

GLOBAL  (SAFE) 

LOCAL  (FEM) 

ill  > 

GLOBAL  (SAFE) 

Figure  8  Schematic  of  global-local  approach. 


The  SAFE  method  is  employed  to  calculate 
the  guided  wave  modeshapes  O  and 
eigenvalues  £,  for  any  infinitely-uniform 
cross-section.  These  modeshapes  and 
eigenvalues  are  used  to  calculate  the  forces 
and  displacements  on  the  boundary 
separating  the  Global  from  the  Local  region. 
In  the  SAFE  method,  at  each  frequency  go  a 
discrete  number  of  guided  modes  is  obtained. 


i  Mono-d  imensio  na! 

Element 


u,. 


u 


u,. 


U, 

< 


A 


u, 


y2 


■U, 


y3 


For  the  given  frequency,  each  mode  is 

characterized  by  a  wavenumber,  §,  and  by  a  F'gure  9  Discretization  in  the  semi-analytical 
displacement  distribution  over  the  cross-  element  method. 

section.  For  a  Cartesian  reference  system,  the  waveguide  cross-section  is  set  in  the  y-z  plane  while 
the  x-axis  is  parallel  to  the  waveguide  length  (Figure  9). 


Subdividing  the  cross-section  via  finite  elements,  the  approximated  displacement  at  a  point, 
u  =  u(x,y,z,t) ,  is  given  as: 


it  =  NUL’(J^X  (17) 

where  N  =  N(y,z)  is  the  matrix  of  the  shape  functions,  Ue  is  the  nodal  displacement  vector  for  the 
e-th  clement,  /  is  the  time  variable  and  /  the  imaginary  unit.  It  can  be  noted  that  the  displacement  is 
described  by  the  product  of  an  approximated  finite  element  field  over  the  waveguide  cross-section 
with  exact  time  harmonic  functions,  ,  in  the  propagation  direction,  x.  The  compatibility  and 

constitutive  equations  can  be  written  in  synthetic  matrix  forms  as: 
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£  =  Du 


a  =  C‘e 


(18) 


where  £  and  G  arc  the  strain  and  stress  vectors  respectively,  D  is  the  compatibility  operator  and 
C  is  the  stiffness  tensor  which  is  general  can  be  complex  (viscoelastic  behavior).  More  details  on 
the  compatibility  operator  can  be  found  in  (Gopalakrishnan  and  Doyle  1995)  and  (Doyle  1995).  The 
principle  of  virtual  work  with  the  compatibility  and  constitutive  laws,  Eqs.(17)  and  (18),  leads  to  the 
following  energy  balance  equation 

J  5uTr  dr  =J  SuT(pii)dV+j  S(uD)T  C’DudV  ( 1 9) 

where  r  is  the  waveguide  cross-section  area,  V  is  the  waveguide  volume,  r  is  the  external  traction 

vector  and  the  overdot  means  time  derivative.  The  finite  element  procedure  reduces  Eq.(19)  to  the 
set  of  algebraic  equations: 


[ASB\U  ®  =  P  (20) 

where  the  subscript  2M  indicates  the  dimension  of  the  problem  with  M  the  number  of  total  degrees 
of  freedom  of  the  cross-sectional  mesh.  Details  on  the  complex  matrices  A  and  B  can  be  found  in 
Hayashi  et  al  2003.  By  setting  p- 0  in  Eq.(20),  the  associated  eigenvalue  problem  can  be  solved  as 
^(co)  which  gives  the  wavenumbers  for  the  different  Lamb  modes  of  the  structure  at  frequency  co. 
For  each  frequency  (;),  2M complex  eigenvalues  tJV,  and  2M complex  eigenvectors  Om  arc  obtained. 
The  solution  is  symmetric,  i.e.  for  each  pair  (<^,  ®m),  representing  a  forward  guided  mode,  a  pair 
exists  representing  the  corresponding  backforward  mode.  The  first  M  components  of  Om  describe 
the  cross-scctional  mode  shapes  of  the  m-th  mode.  Once  ^  is  known,  the  dispersion  curves  can  be 
easily  computed.  The  phase  velocity  can  be  evaluated  by  the  expression c  h  =  co!  t;reu, ,  where  %real  is 

the  real  part  of  the  wavenumber.  The  imaginary  part  of  the  wavenumber  is  the  attenuation, 
att=  E,jmag ,  in  Nepers  per  meter. 

3.2.  Global-Local  Method 

The  global  functions,  represented  by  <Dn,  arc  the  displacements  that  satisfy  the  governing  equations 
in  the  absence  of  the  discontinuity.  The  displacements  at  the  mesh  boundary  are  written  as  a  sum  of 
the  contribution  due  to  the  incident  wave  and  that  due  to  the  scattered  waves.  This  displacement  is 
then  further  used  to  satisfy  the  continuity  conditions  at  the  boundary. 

The  displacements  at  the  left  boundary  are  written  as: 


K incident  ^ scattered 


7=1 


(21a) 
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»,  =  ^®;«*'&")+[c'][D']  (2 lb) 

where  the  time-harmonic  dependence  is  implicit. 

Rattanawangeharoen  et  al  (1997)  gives  a  detailed  description  of  matrices  [G  •  In  essence, 

rG”J  is  a  matrix  consisting  of  the  modeshapes  of  the  scattered  modes  whereas  [D"J  is  a 

column  vector  which  contains  the  resultant  displacement  due  to  these  modes  on  the  left 
boundary. 

In  equations  (21a,  21b),  d>  represents  the  modeshape  of  the  incident  waves  travelling  to  the  right, 
<!>'  represents  the  modeshapes  of  the  reflected  waves  travelling  to  the  left;  Am ,  O) ,  and 

represent  the  amplitude,  modeshape  and  wavenumber  of  the  incident  wave  respectively;  A~ 

represents  the  unknown  amplitude  of  the  reflected  waves;  superscript  *+’  represents  a  wave 
travelling  in  the  right  and  a  superscript  represents  a  wave  travelling  in  the  left  direction;  x' 
represents  the  absolute  value  of  the  distance  of  the  left  boundary  from  the  center  of  the  meshed 
region;  L  is  the  total  number  of  modes  considered.  In  the  same  way,  the  displacements  at  the  right 
boundary  are  written  as: 

^ incident  '  ^ scattered 


(22a) 

(22b) 

In  equations  (22a)-(22b),  A]  represents  the  unknown  amplitude  of  the  transmitted  waves  and  x 

represents  the  absolute  value  of  the  distance  of  the  right  boundary  from  the  center  of  the  meshed 
region. 

The  forces  at  the  boundaries  can  be  calculated  by  computing  the  consistent  nodal  force  vectors  due 
to  displacements  (21b)  and  (22b).  The  nodal  forces  at  the  left  boundary  arc: 

/=-A,4V<-«*  '+[f-][d  ]  (23) 

Fj  is  the  consistent  nodal  force  vector  calculated  from  the  displacement  veetor  O,.  The  nodal  forces 
at  the  right  boundary  are: 

X=-4^>'<°')+[/r+][^+]  (24) 


Distribution  B:  Approved  for  public  release;  distribution  is  unlimited. 


13 


3.3.  Finite  Element  Model  for  Local  Region 

The  conventional  discretization  process  in  the  finite  element  method  yields: 


*{j'XsK«HfcX,,.}=o 


(25a) 


Mr  =<fc}rfe»}r) 


(25b) 


(25c) 


Here,  [AT],  and  [M]  are  global  stiffness  and  mass  matrices  respectively;  co  is  the  haimonie 

frequency;  {g,}  and  {gfi}  are  the  nodal  displacement  vectors  corresponding  to  interior  and 

boundary  nodes  respectively;  {P„}  is  the  interaction  force  vector  at  the  boundary  nodes;  and  {</ } 

denotes  complex  conjugate  of  {<y}.  Forces  on  the  boundary  nodes  due  to  the  entire  FE  region,  when 

displacement  is  imposed  only  on  the  boundaries,  can  be  calculated  by  performing  static 
condensation  on  [5]  such  that: 


[s;,]=[5„]-[5„][s„r'[s„]  (26) 

3.4.  Global  Solution 

The  nodal  forces  on  the  boundary  due  to  displacements  (21b)  and  (21b),  as  contributed  from  the  FE 
portion,  are: 


For  traction  continuity,  from  equations  (21),  (22),  and  (27): 


I 


{/.}} 

{fr}\ 


rFES4 
J  boundary 


(28) 


Rearranging  the  terms,  the  following  final  equation  is  obtained: 
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(*;  Hs»  ]  {«*  i 


WJ- 


A.Ke,<ri-n 

a.  f:Al'  * 


(29a) 

(29b) 

(29c) 


{D}r=({D-/{D*}r) 


(29d) 


[F-] 

o  [r] 


M 

0  [O’] 


(290 


Equation  (29a)  is  solved  using  least  squares  to  calculate  the  coefficients  A~ ,  and  A * .  Reflection 
coefficients  (/?„)  and  transmission  coefficients  (71)  (due  to  the  ith  incident  mode)  are  defined 


below: 


(30) 

(31) 


3.5  Energy  Conserv  ation 

Energies  are  carried  only  by  the  propagating  modes,  which  include  incident,  reflected  and 
transmitted  modes.  The  time  averaged  values  of  energy  flux  associated  with  the  j*  reflected  and 
transmitted  modes  through  the  cross-section  is  given  by: 
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where 


(33) 


is  the  participation  factor  for  the  j"  propagating  mode  and  the  integration  is  carried  over  the 
thickness  of  the  plate  ( h ). 

3.6.  Numerical  Results  And  Experimental  Comparisons 

3.6.1.  Reflection  for  Notched  Aluminum  Plate 


This  study  examined  a  semi-infinite  plate  with  a  notch  defect.  The  plate  modeled  was  aluminum, 
1.58  mm  (1/16  in)  in  thickness.  The  first  notch  considered  was  0.16  mm  in  width  and  0.32  mm  in 
depth  (Figure  10(top)).  A  finite  region  surrounding  the  notch  (local  domain)  was  discretized  by 
rectangular,  four-node  finite  elements.  A  total  of  ten  elements  were  used  in  the  thiekness  of  the 
plate,  and  the  mesh  was  kept  uniform  throughout  the  local  domain.  SAFE  was  used  for  the  global 
domain.  The  zero-order  symmetric  mode,  So,  was  considered  as  the  incident  mode. 

The  results  are  shown  in  the  plot  of  Figure  10(top)  in  terms  of  reflected  energy  relative  to  the 
incident,  unit  energy.  As  expected,  the  notch  causes  So-Ao  mode  conversion  since  it  breaks  the 
symmetry  of  the  waveguide.  As  a  check  for  the  accuracy  of  the  analysis,  energy  conservation  is 
preserved  as  the  reflected  So  and  Ao  energies  sum  to  the  incident  So  energy.  Mode  conversion  is 

found  more  pronounced  at 
particular  frequency  values, 
namely  500  kHz,  730  kHz,  880 
kHz  and  950  kHz. 

A  more  shallow  discontinuity 
was  examined  next.  This  was  a 
notch  extending  1 .6  mm  in  width 
and  0.16  mm  in  depth,  which 
more  closely  represented 
corrosion  damage  (Figure 
10(bottom)).  It  can  be  seen  that 
the  resulting  reflection  spectra  are 
different  for  the  two  defect  cases. 
This  suggests  the  potential  for 

quantitative,  rather  than 
Figure  1 0  Reflection  spectra  for  notch  defects  in  A1  plate.  qualitative  defect  detection. 
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3.6.2.  Reflection  and  Transmission  for  Notched  Aluminum  Plate 


Figure  11  Reflection  and 
transmission  spectra  for  notch 
defects  in  an  A1  plate. 

In  this  case  the  global-local  method 
was  used  to  predict  the  scattered  field 
in  both  reflection  and  transmission. 
The  plate  modeled  was  aluminum, 
1.6  mm  in  thickness  with  a  half 
thickness  notch  in  the  middle. 
Transmitted  and  reflected  So  and  Ao 
modes  were  computed  for  an  incident 
So  mode.  Energy  conservation  was 
again  checked  to  test  the  validity  of 
the  method  (Figure  11). 


3.6.3.  Time  Histories 


The  global-local  method  was 
also  used  to  predict  the  time- 
domain  waveforms  in  an 
aluminum  plate  subjected  to 
broadband,  pulsed  laser 
excitation.  The  response  was 
experimentally  measured  for 
the  case  of  a  plate  1 .58  mm  in 
thickness  with  a  sensor- 
excitation  distance  of  20  cm. 

Figure  12(a)  shows  the 

geometry  of  the  problem. 

Figure  12(b)  and  Figure  12(c) 
show  the  comparison  of  the 
theoretical  predictions  with  experimentally  measured  results.  The  numerical  predictions  were 
calculated  by  inverting  the  frequency  response  calculated  from  the  SAFE-based  global-local 
method. 

It  can  be  seen  that  a  small  discrepancy  exists  between  the  model  and  the  experiment  in  the  wave 
arrival  time  of  So,  whereas  good  match  is  found  for  Ao.  The  So  discrepancy  may  be  due  to  an 
imperfect  knowledge  of  the  material  in  the  test.  It  is  nevertheless  worth  noting  that  both  the 
numerical  predictions  and  the  experimental  results  show  a  decrease  in  the  amplitude  of  the  Ao  mode 
in  the  case  of  the  notched  plate.  The  numerical  prediction  also  shows  significantly  more  mode 
conversion  when  compared  to  the  experimental  results.  This  is  primarily  because  of  energy 
conservation  concerns  for  the  plain  strain  approximation  of  a  three  dimensional  problem.  While 
energy  is  allowed  to  scatter  in  random  directions  in  the  experimental  case,  reducing  the  problem  to  2 
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Figure  12  Time  histories  for  notched  Al  plate:  experiments 
and  simulations. 
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dimensions  forees  the  scattered  energy  from  the  Ao  mode  to  be  transferred  to  the  So  mode.  Further 
studies,  including  adding  higher-order  scattered  modes,  refining  the  finite  clement  mesh,  and 
extending  the  2-D  global-local  formulation  to  3-D,  are  being  conducted  to  obtain  a  better 
convergence  of  the  model  with  the  experiments. 

4.  INSONIFICATION  APPROACH 

4.1.  Chaotically-modulated  Insonification  and  Time  Series  Predictive  Modeling 

These  GUW  techniques  are  well  established  and  can  work  in  certain  structures  with  simple 
geometries  (e.g.,  plates,  beams,  ete.)  or  on  sections  with  constant  cross-seetion  properties  in  the 
wave  propagation  direction  (e.g.,  rails).  However,  these  methods  cannot  as  easily  be  applied  to 
irregular  geometries,  such  as  bolted  joints,  stiffened  ribs,  or  curve  panels,  because  of  mode 
conversion  and  wave  interference  effects  that  arise  as  a  result  of  the  complex  interfaces  and 
mechanical  impedanee  mismatches  in  these  constructions.  Instead,  some  researchers  have 
attempted  to  employ  bulk  insonification,  where  an  ultrasonic  source  is  excited  and  the  resultant 
long-time,  or  diffuse,  wave  field  is  examined  to  identify  structural  changes  (Michaels  and 
Michaels,  2005).  This  method  is  complementary  to  GUW  mode-locked  analyses  for  structures 
w'ith  complex  boundary  conditions  or  geometries  that  make  tracking  and  analysis  of  a  single 
propagating  mode  difficult  or  impossible. 

The  new  approach  created  here  uses  a  ehaotieally-modulated  ultrasonic  wave  as  the  probe. 
Methods  that  employ  chaotic  excitations  and  attractor-based  prediction  error  algorithms  (but  not 
in  the  ultrasonic  domain)  have  demonstrated  the  eapaeity  to  detect  defeets  in  various  test  bed 
structures  with  enhanced  sensitivity  over  traditional  vibration-domain  analyses  (Nichols,  Todd 
and  Wait,  2003;  Todd,  Erickson,  Chang,  Lee  and  Niehols,  2004).  The  key  lies  in  the  phase-space 
diversity  of  the  chaotic  sideband  structure.  Dispersion  issues  are  less  important,  as  the  primary 
point  is  to  use  the  phase-diverse  (chaotic)  waveform  as  a  constant  probe  souree  and  employ  some 
form  of  waveform  comparison  or  time  series  predictive  model  to  contrast  a  “baseline”  condition 
with  a  test  condition.  The  ehaotie  ultrasonic  waves  are  fundamentally  created  via  amplitude 
modulation,  i.c.,  by  multiplying  a  single  ultrasonic  frequency  tone  by  an  amplitude  envelope  that 
is  created  by  a  separate  chaotic  process.  A  chaotic  waveform  is  able  to  enhance  prediction  error 
based  features  beeause  the  signal  is  deterministic;  it  is  not  just  frequency  structure  that’s 
important.  Using  a  phase-randomized  signal  with  frequency  content  equivalent  to  a  chaotic 
waveform  results  in  extracted  features  that  are  not  able  to  identify  damage  as  well  as 
deterministic  chaotic  waveforms  (Niehols,  Todd  and  Wait,  2003).  The  waveform  appears  as  a 
narrowband,  chaotically-modulated  signal  centered  at  some  central  carrier  frequency.  Although 
not  required,  it  is  usually  best  to  select  a  carrier  frequency  where  dispersion  is  still  somewhat 
minimal  relatively  minor  dispersion  that  oeeurs  in  aluminum  at  that  frequency  (although  this  is 
not  a  requirement),  beeause  of  the  bar  thickness  that  is  used  in  the  experimental  study,  and 
beeause  the  traveling  modes  are  well-separated  in  phase  velocity  space. 

Creation  of  this  signal  is  a  multi-step  proecss  and  starts  with  the  generation  of  a  1  Hz  sine  wave 
with  a  time  step  dt=l/fs,  where  fs  is  the  sampling  frequency  (chosen  to  be  4  MHz  for 
experimental  investigation  below  due  to  data  acquisition  limitations).  A  1  Hz  sine  wave  is  used 
because  a  standard  Lorenz  chaotic  signal  has  significant  frequency  information  only  in  the  1  Hz 
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range  and  therefore  when  the  ehaotie  signal  is  used  to  window'  the  sine  wave  they  have  similar 
frequency  eontent.  At  the  end  of  the  w  indowing  proeess  the  new  windowed  chaotic  signal  ean 
then  be  shifted  to  the  desired  center  frequency  by  changing  the  size  of  the  time  step  used  in  the 
time  vector.  The  chaotic  signal  is  constructed  using  the  output  of  the  x  variable  from  the 
following  three-dimensional  nonlinear  Lorenz  system: 

x  =  10(y- jr) 
y  =  (—xz  +  28  x-y) 

z=(xy-Sz/3)  (34) 

There  is  nothing  unique  about  the  Lorenz  system  for  generating  chaotic  output;  any  system 
capable  of  producing  a  chaotic  output  is  suitable,  but  the  Lorenz  system  shown  in  (34)  has  a 
robust  parameter  region  for  producing  ehaotie  output  and  was  selcetcd  for  this  study.  Equation 
(34)  is  integrated  using  a  time-step  dt*R ,  where  R  is  a  frequency  ratio  that  can  be  modified  to 
change  the  frequency  bandwidth  of  the  ehaotie  signal.  For  this  study  we  use  several  values  of  R 
which  affect  the  frequency  regime  in  w'hieh  the  power  of  the  chaotic  signal  lies.  For  example,  a 
value  of  R= 1/3  ereates  a  signal  in  which  the  significant  frequency  information  (as  determined  by 
a  loss  of  40  dB  in  power  spectral  density)  is  less  than  1.5  Hz.  A  value  of  /?=1  /30  would  result  in 
significant  frequency  information  being  less  than  0.2  Hz.  This  ehaotie  signal  is  normalized 
through  division  by  the  maximum  of  the  absolute  value  of  the  signal  so  that  the  values  range 
from  -1  to  l.A  modulated  signal  is  then  created  using  the  following  equation: 

y  =  sin(2;r*  t)*  (1  +  d*  x)  ^5) 

where  x  is  the  Lorenz  waveform,  sin(27r*f)  is  the  originally  created  sine  wave,  d  is  the 
modulation  depth,  and  y  is  the  chaotically  amplitude-modulated  output  signal.  If  the  value  of  d, 
whieh  controls  signal  bandwidth,  is  greater  than  one,  the  resulting  signal  will  be  over-modulated 
and  will  result  in  a  phase  inversion  at  the  points  where  |c/*x|>l.  These  phase  inversions  would  be 
detrimental  to  any  prediction  algorithm,  and  d  is  therefore  restricted  to  the  range  0<J<1.  The 
amplitude  modulated  signal  is  then  upconvertcd  to  the  target  earner  wave  frequency  by 
multiplying  dt*fc.  The  resulting  ehaotie  time  series  are  also  smoothed  in  time  at  its  boundaries 
with  a  trapezoidal  window  to  facilitate  launching  with  piezoelectric  devices.  The  CUWs  were 
launched  over  a  2  ms  time  period.  Figure  1 3  depicts  the  actual  time  series  for  various  parameter 
combinations,  and  Figure  14  shows  the  effect  of  changing  the  frequency  ratio  R  as  well  as  the 
modulation  depth  d  on  the  power  spectral  density  of  the  modulated  carrier  wave.  The 
effectiveness  of  the  damage  detection  technique  used  in  this  study  is  highly  dependent  on  carrier 
frequency  and  frequency  ratio  but  appears  to  be  relatively  insensitive  to  modulation  depth, 
provided  a  value  significantly  larger  than  zero  is  used  (Fasel,  Olson,  and  Todd,  2008). 
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Figure  13  Modulated  signals  using  (top  left)  R  =  0.10  and  d  =  0.5;  (top  right)  R  = 
0.10  and  d=  1.0;  (bottom  left)  R  =  0.33  and  d  =  0.5;  (bottom  right)  R  =  0.33  and  d 
=  1.0. 


30 


Figure  14  (left)  Power  spectral  densities  of  modulated  signals  using  d=0.5  and 
(right)  power  spectral  densities  of  modulated  signals  using  R=0.33. 

Once  the  CUW  has  been  created  it  is  then  imparted  to  the  structure  through  an  MFC  (or  similar 
piezoelectric)  sensor/actuator,  the  characteristics  of  which  were  described  in  Section  2.  The 
waveform  is  detected  by  a  second  sensing  MFC  in  a  pitch/catch  mode.  The  primary  task  at  this 
point  in  the  damage  detection  scheme  is  to  decide  what  feature(s)  from  the  measured  waveform 
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may  be  extracted  to  best  assess  the  condition  status  of  the  structure.  A  novel  statistical 
classification  technique  with  its  basis  in  information  theory  is  employed  for  this  study. 

A  fundamental  theorem  of  Shannon's  information  theory  says,  paraphrased  intuitively,  “the  best 
compression  for  any  given  data  set  comes  from  a  codebook  designed  exactly  for  the  statistics  of 
that  souree,  any  other  eodebook  will  give  worse  results”  (Shannon  and  Weaver,  1949).  For 
instance,  if  one  has  a  eodebook  (e.g.,  taking  language  elements  like  words  into  shorter  eodcs) 
consisting  of  English  words  and  another  consisting  of  French  words,  then  a  new  time  series  of 
letters  ean  be  represented  in  the  shortest  compressed  format  when  using  the  English  eodebook 
versus  all  others  if  the  new  text  is,  in  faet,  written  in  English.  Compression  performance  is  the 
elassie  text  categorization  methodology.  Modem  statistics  has  melded  ideas  of  information 
theory  to  extend  to  continuous  signals,  where  compression  performance  is  intimately  tied  with 
out-of-sample  prediction  error,  and  a  eodebook  is  the  model  for  a  source  (something  that 
produces  time  series  data).  In  this  ease  the  souree  arises  from  an  actual  physical  proeess  (the 
guided  wave  propagating  through  a  bolted  joint).  This  idea  leads  to  a  procedure  for  classifying 
time  series  using  eross-prediction  error,  as  literal  “data  compression”  is  not  actually  neeessary, 
just  its  “virtual”  performance.  First,  a  suitable  model  of  the  measured  response  time  series  is 
required,  as  opposed  to  the  eomputationally-intensive,  first-principles  physical  spatially-extended 
finite  element  model  of  the  joint.  The  perfect  model  for  classification  is  not  needed,  as  it  can 
work  well  with  reasonable  model  misspeeifieation.  However,  the  better  the  underlying  statistical 
model  is,  of  course,  the  more  the  classification  performance  will  improve. 

The  first  step  in  this  method  employs  the  use  of  autoregressive  (AR)  models,  which  have 
previously  been  shown  to  be  useful  in  damage  classification  schemes  (Sohn  and  Farrar,  2001; 
Sohn,  Farrar.  Hunter  and  Worden,  2001).  The  discretely  observed  output  time  series  x(n)  is 
modeled  with  an  AR  model  of  the  form 


p 

\(n)  =  Y,0CjX(n  -  /)  +  e(/i) 

(36) 

where  p  is  the  order  of  the  AR  model  with  associated  eoeffieients  a,  and  residual  error  e(«).  In 
this  study  it  was  determined  that  an  order  of  p= 25  effectively  models  the  sensed  waveforms.  The 
AR  coefficients  are  estimated  through  minimization  of  the  sum  of  squared  forward  prediction 
errors  (Broekwell  and  Davis,  1991).  All  signals  are  normalized  by  dividing  the  standard 
deviation  of  the  signal  before  use  of  the  AR  model. 

The  proeess  behind  the  classification  technique  is  as  follows.  First,  a  set  of  distinct  input  signals 
are  ereated  from  the  data-generating  proeess  that  has  been  previously  described.  For  eaeh  of 
these  input  signals  a  structural  response  is  recorded  under  various  damage  conditions  of  interest. 
For  eaeh  of  these  responses  AR  coefficients  are  estimated  using  the  above  outlined  method.  The 
sets  of  AR  eoeffieients  for  eaeh  damage  condition  and  eaeh  input  signal  forms  a  training 
database.  Eaeh  configuration  from  a  test  signal,  observed  when  the  system  was  driven  by  an 
input  that  was  not  ever  observed  in  the  database  before  (yet  eame  from  the  same  data-generating 
proeess  as  in  the  original  information  theoretic  context)  is  classified  using  these  coefficients.  A 
new  input  signal  (created  from  the  same  underlying  proeess  as  the  database  input  signals)  is  then 
applied  to  the  structure  when  the  structure  is  in  an  unknown  damage  state.  Eaeh  set  of  AR 


Distribution  B:  Approved  for  publie  release;  distribution  is  unlimited. 


21 


Fit  \R 
models  to 
each 
response 


Interrogate 
structure  in 
unknown 
tnage  slat 


Save  A  R 
coefficient 
“training” 
database 


Analyze  test 
signals  and 
predict 
damage  state 


coefficients  in  the  training  database 
for  the  first  of  the  input  signals  (that 
is  each  set  of  AR  coefficients  which 
describe  different  damage  conditions 
for  the  first  input  signal)  is  then  used 
to  estimate  the  structural  response  to 
the  new  input  signal.  One  set  of 
coefficients  from  the  training 
database  will  minimize  the  sum  of 
the  squared  residual  errors.  The 
structural  condition  state  that  the 
structure  was  in  when  these  AR 
coefficients  were  recorded  is  scored 
as  the  “vote”  for  the  unknown 


Figure  15  Classification  paradigm. 


condition  using  that  particular  input 
signal.  This  comparison  then  takes 
place  for  each  of  the  remaining  input  signals  in  the  training  database.  This  entire  process  is  then 
repeated  using  multiple  input  signals  that  are  imparted  to  the  structure  in  its  unknown  state.  The 
votes  for  each  condition  are  then  summed  and  the  condition  with  the  plurality  of  votes  is  the 
estimated  condition  of  the  structure.  The  statistical  classification  paradigm  can  be  shown  visually 
in  Figure  15. 


4.2.  Application  to  Bolt  Preload  Loss  Detection  in  Aluminum  Structural  Joints 


Two  experimental  metallic  test  platforms  employing  fastened  (bolted)  joints  were  created  to  test 
the  effectiveness  of  this  method  using  real  materials  and  the  accompanying  experimental  noise. 
This  noise  level  is  minimized  to  the  greatest  extent  possible  by  employing  shielded  cabling  as 
well  as  moving  the  test  structures  an  acceptable  distance  awray  from  the  signal  amplifier,  which 
emits  electromagnetic  interference  (EMI)  during  the  actuation  proeess.  Each  input  signal  is  also 
applied  to  the  structure  50  times  and  then  averaged  in  a  further  attempt  to  reduce  experimental 
noise.  The  actuation  signal  is  created  by  the  output  channel  of  a  National  Instruments  PCI-61 10 
DAQ  card  at  a  rate  of  4  MHz  and  routed  through  a  Krohn-Hite  7602  wideband  power  amplifier. 
This  amplified  signal  is  sent  to  the  actuation  MFC  while  the  sensing  MFC  simultaneously 
samples  the  structural  response  at  a  rate  of  4  MHz. 

Single-bolt  lap  joint 

The  first  experimental  apparatus  on  which  testing  was  carried  out  is  the  single  bolt  lap  joint 
shown  in  Figure  16.  The  structure  is  made  up  of  two  aluminum  bars  (0.3m  x  .05m  x  9.5mm) 
connected  to  eaeh  other  with  a  single  bolt.  Two  MFC  patches  were  attached  to  the  structure  with 
one  on  both  sides  of  the  joint.  Each  of  these  Smart  Material  Corp.  MFC  patches  (M  2814  P2) 
have  an  active  area  of  28mm  x  14mm,  are  approximately  0.3mm  thick  and  are  bonded  to  the 
structure  50mm  from  the  bolted  connection  (on  eaeh  side)  using  cyanoacrylate.  Due  to  the 
symmetry  of  the  problem  it  is  not  important  w'hieh  patch  is  used  as  the  sensor  and  which  as  the 
actuator  and  either  configuration  will  yield  similar  results.  For  this  experiment  the  transducer  on 
the  left  side  of  the  lap  joint  in  Figure  16  was  used  as  the  actuator. 
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Figure  16  Single-bolt  lap  joint  test  article. 


Table 


:  Classification  "vote"  distribution  of  experimental  lap  joint  data 
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24 

122 

79 
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In  this  study,  data  were  taken  at  each  step  of  a  bolt  tightening  sequence  in  which  the  bolt 
condition  is:  'loose'  (Condition  1),  'finger-tight'  (Condition  2),  3.5  N-m  (Condition  3),  and  14  N- 
m  (Condition  4).  This  sequence  is  then  repeated  three  times  to  simulate  assembly  and 
disassembly  of  the  joint  in  a  real  structure.  The  first  two  assembly/disassembly  sequences  were 
used  to  create  a  training  database.  Structural  responses  obtained  during  the  third  sequence  were 
used  as  test  inputs.  Table  1  shows  the  vote  results  for  the  4  conditions.  Again,  if  all  individual 
test  inputs  are  classified  correctly,  votes  would  appear  only  along  the  diagonal  in  bold. 

The  true  bolt  condition  was  correctly  assessed  by  the  statistical  classification  algorithm  in  all 
cases  but  Condition  4  (the  most  tight),  which  was  estimated  by  vote-counting  to  be  in  Condition 
3.  There  arc  several  factors  contributing  to  this  damage  case  misidentification.  First,  specifying 
bolt  torque  on  a  real  joint  can  be  difficult  and  in  this  experiment  an  inexpensive  torque  wrench 
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with  a  fairly  low  resolution  was  used.  Second,  the  transfer  relationship  between  torque  and 
preload  is  hysteretic,  nonlinear  and  is  highly  dependent  on  local  contact,  which  will  vary  each 
time  the  bolt  is  tightened.  Third,  in  this  experiment  it  was  difficult  to  maintain  the  boundary 
conditions  of  the  lap  joint  between  tests  and  it  is  believed  that  this  also  led  to  inflated  number  of 
incorrect  votes.  Given  these  concerns,  the  actual  preload  indicated  by  a  particular  torque  level 
may  vary  significantly  from  test  to  test  and  almost  certainly  contributed  to  a  mueh  lower 
percentage  of  correct  identification  of  individual  test  cases  (66%)  than  were  seen  in  simulation 
results  performed  separately. 

In  future  tests  an  instrumented  bolt  will  be  used  so  that  a  direct  measure  of  preload  will  be 
available  instead  of  just  bolt  torque.  This  should  improve  results,  however  in  real  world 
situations  bolt  preload  will  be  specified  by  torque  specifications.  Other  improvements  that  are 
not  dependent  on  knowing  exact  bolt  preload  level  are  possible.  Foremost  among  these  arc  the 
choice  of  parameters  that  affect  the  creation  of  the  input  time  signals  (carrier  frequency  fc, 
frequency  ratio  R,  modulation  depth  d)  as  well  as  feature  extraction  (AR  model  order,  size  of 
training  and  test  databases).  Using  genetic  algorithms  (specifically  differential  evolution)  to 
create  an  optimal  input  waveform  for  maximum  damage  discernment  has  already  been 
investigated  (Fasel,  Olson  and  Todd,  2008).  This  method  showed  significant  improvement  (two 
orders  of  magnitude)  in  solution  ‘fitness’  over  a  random  set  of  input  parameters  (which  is  what 
was  used  in  this  study).  This  statistical  classification  method  has  also  been  used  on  a  composite 
plate-to-spar  bond  with  multiple  disbond  sizes  as  well  as  a  poorly  cured  bond.  Results  from  these 
experiments  have  not  yet  been  published  but  show  that  the  ability  of  this  method  to  eorreetly 
identify  damage  state  is  highly  dependent  on  carrier  frequency  fc  and  frequency  ratio  R.  As  well, 
this  classification  scheme  might  be  best  used,  in  the  ease  of  detecting  bolt  preload  loss,  by 
specifying  a  critical  threshold  value  of  preload  level  above  which  the  joint  is  considered  healthy 
and  below  which  the  joint  is  considered  damaged.  This  critical  threshold  value  is  further 
examined  in  the  second  experimental  structure  being  investigated  in  this  study. 

Multiple-bolt  portal  structure 

As  mentioned  in  the  previous  subseetion,  it  is  believed  that  a  test  bed  with  more  reliable  end 
boundary  conditions  should  result  in  a  greater  percentage  of  correet  classifications.  It  was  also 
desired  to  test  a  structure  that  had  multiple  bolted  connections  in  order  to  examine  the  ability  of 
the  statistical  classification  algorithm  to  identify  multiple  damage  locations  within  a  structure. 
This  ability  to  locate  damage  within  a  multiple  bolt  structure  is  a  necessity,  as  virtually  all 
practical  field  applications  will  fall  into  this  category.  Therefore,  the  aluminum  frame  structure 
shown  in  Figure  17  was  employed  to  address  these  concerns.  The  two  side  bars  are  the  same 
dimensions  as  the  bars  used  in  the  single  bolt  lap  joint  experiment  (0.3m  x  .05m  x  9.5mm)  and 
the  top  bar  is  twice  as  long  (0.6m  x  ,05m  x  9.5mm).  The  actuating  MFC  is  0.08m  from  the  close 
end  of  the  bar  (next  to  bolt  1)  and  0.45m  from  the  far  end  of  the  bar  (next  to  bolt  2).  The 
actuating  MFC  is  placed  asymmetrically  to  remove  the  symmetry  of  the  structure  and  to  make 
individual  bolt  damage  state  identification  more  possible.  Sensing  MFC  1  and  MFC  2  are  0.08m 
from  the  top  of  their  respective  side  bars.  The  bolts  are  each  0.4m  from  the  end  of  bar. 
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Figure  17  Multiple-bolt  test  structure. 


Table  2  shows  the  damage  cases  that  were  considered  in  this  study.  'Tight'  indicates  120  in-lb, 
'finger  tight'  indicates  nominal  preload  (less  than  30  in-lb),  and  'loose'  indicates  no  preload. 
While  there  are  thus  7  “conditions”  defined,  as  indicated,  there  are  only  3  damage  levels  for  each 
bolt.  The  last  7  eases  were  used  as  test  cases  against  the  training  database  created  using  the  first 
14  eases. 

Table  2:  Test  conditions  of  the  multiple-joint  structure 


Case 

Bolt  1  Condition 

Bolt  2  Condition 

1 

Tight 

Tight 

2 

Finger  Tight 

Tight 

3 

Loose 

Tight 

4 

Tight 

Finger  Tight 

5 

Tight 

Loose 

6 

Finger  Tight 

Finger  Tight 

7 

Loose 

Loose 

8 

Tight 

Tight 

9 

Finger  Tight 

Tight 

10 

Loose 

Tight 

11 

Tight 

Finger  Tight 

12 

Tight 

Loose 

13 

Finger  Tight 

Finger  Tight 

14 

Loose 

Loose 

15 

Tight 

Tight 

16 

Finger  Tight 

Tight 
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17 

Loose 

Tight 

18 

Tight 

Finger  Tight 

19 

Tight 

Loose 

20 

Finger  Tight 

Finger  Tight 

21 

Loose 

Loose 

Table  3:  Classification  "vote"  distribution  of  multiple 


■joint  frame  data 


Damage 

Case 

MFC  l(Bolt  1) 

MFC  2  (Bolt  2) 

Tight 

Finger 

Tight 

Loose 

Tight 

Finger 

Tight 

Loose 

15 

225 

0 

0 

225 

0 

0 

16 

0 

152 

73 

225 

0 

0 

17 

2 

115 

108 

225 

0 

0 

18 

225 

0 

0 

0 

134 

91 

19 

225 

0 

0 

0 

98 

127 

20 

1 

140 

84 

0 

196 

29 

21 

0 

3 

222 

06 

1 

224 

The  vote  chart  for  each  MFC  sensor  is  shown  in  Table  3.  The  bold  numbers  in  each  row  indicate 
the  true  condition  of  the  bolt.  Therefore  a  correct  classification  is  made  if  the  bold  number  is  the 
largest  in  its  row.  As  such,  the  correct  classification  was  made  in  each  case  except  for  bolt  1  in 
damage  case  17.  The  damage  localization  ability  of  this  method  is  still  strong  as  the  overall 
percentage  of  correctly  identified  individual  test  cases  is  84%.  The  'tight'  condition  was  classified 
correctly  for  almost  every  individual  test  signal.  However,  the  distinction  between  'finger  tight' 
and  'loose'  is  less  clear.  This  unclear  discemability  between  the  ‘finger  tight’  and  ‘loose’  damage 
conditions  invites  the  employment  of  a  critical  threshold  value  as  discussed  previously. 


Table  4:  Classification  Distribution  Multi-joint  Experimental  Data 


Damage 

Case 

MFC  l(Bolt  1) 

MFC  2  (Bolt  2) 

Tight 

Loose 

Tight 

Loose 

15 

225 

0 

225 

0 

16 

0 

225 

225 

0 

17 

2 

223 

225 

0 

18 

225 

0 

0 

225 

19 

225 

0 

0 

225 

20 

1 

224 

0 

225 

21 

0 

225 

0 

225 

Therefore,  the  categories  'finger  tight'  and  'loose'  were  combined  into  a  more  general  'loose' 
category  by  establishing  the  critical  threshold  value  at  a  preload  level  of  ‘finger  tight’.  In  this 
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attempt  to  make  a  purely  healthy/unhealthy  joint  status  determination,  proper  classification  is 
achieved  with  each  damage  ease.  This  simple  classification  works  so  well  that  votes  for 
individual  test  responses  ehoose  the  eorreet  joint  configuration  greater  than  99%  of  the  time,  as 
ean  be  seen  in  Table  4. 

4.3.  Chaotically-modulated  Insonification  with  Embedology  and  State-space  Modeling 

A  generalization  of  the  time  series  AR  predictive  modeling  approach  in  4.1  is  to  reconstruct 
(from  measured  data)  representations  of  the  state  space  and  use  the  resulting  attractor  geometry 
in  a  predietivc  modeling  eapaeity.  In  praetiee,  it  is  often  difficult  to  measure  experimentally 
each  of  the  state  variables  needed  for  attraetor  reconstruction  (in  faet,  in  a  worst  case,  a  single 
measurement  of  system  response  may  have  been  obtained).  Takens  (Takens  1981)  proved  that, 
under  certain  conditions,  a  state-space  reconstruction  from  a  single  coordinate  is  an  embedding 
of  the  true  state  space  of  the  entire  system.  This  is  accomplished  by  means  of  the  delay- 
coordinate  approach,  whereby  properly  delayed  copies  of  a  single,  discretely-measured 
coordinate  are  used  to  create  an  A/-dimensional  attractor  which  preserves  dynamical  invariants 
(descriptive  quantities  that  relate  to  how  points  are  distributed  on  an  attraetor).  Henee,  the 
attractor  at  time  n  can  be  expressed  as  x„  =  [x(«)x(«  +  r)...,x(«  +  (M  -  l)r)]  where  z  is  an 
appropriate  time  delay.  The  ehoice  of  the  parameters,  Z  and  M,  are  critical  components  to 
properly  embedding  the  attractor.  If  the  delay  is  too  small,  the  resulting  dimension  vectors  will 
be  very  similar  and  will  carry  a  surplus  of  redundant  information.  However,  if  ris  too  large,  the 
coordinates  will  beeome  essentially  unrelated.  Correspondingly,  the  embedding  dimension  must 
be  large  enough  sueh  that  the  reconstructed  orbit  does  not  intersect  itself.  If  the  dimension  is  too 
large,  there  will  be  unnecessary  computational  effort  expended  attempting  to  unfold  an  attraetor 
in  greater  than  M  dimensions  when  only  M  dimensions  are  required.  It  is  an  open  research 
question  about  what  the  best  parameters  to  ehoose  arc,  but  this  work  makes  use  of  the  two  most 
widely  acknowledged  methods  for  choosing  the  proper  delay  in  attractor  reconstruction.  The  first 
method  (Abarbancl  1996)  employs  the  autocorrelation  function  estimate  to  determine  the  delay 
at  which  a  vector  is  least  correlated  with  itself.  Thus,  a  frequently  chosen  delay,  z,  is  taken  to  be 
at  the  first  zero  crossing  of  the  autocorrelation  function.  Another  eommonly  utilized  technique 
for  determining  the  optimal  dimension  is  the  empirieal  false  nearest  neighbors  approach  (Kennel 
and  Abarbancl  1992).  In  this  technique,  the  attractor  is  first  embedded  in  M  =  1  dimensions.  The 
nearest  neighbors  to  each  point  in  a  Euclidean  sense  are  determined.  Next,  the  attraetor  is 
successively  embedded  in  higher  M  values  using  a  previously  chosen  delay,  and  the  number  of 
false  projections  of  nearest  neighbors  is  uncovered.  The  attractor  is  considered  fully  unfolded 
when  the  number  of  false  projections  reaches  zero.  It  should  be  stated  that,  however,  that  when 
embedding  attractors  in  state  spaee  for  SHM  feature  extraction  purposes,  the  goal  is  not 
necessarily  to  best  reconstruct  state-space  geometries.  Rather,  the  objective  is  to  utilize 
reconstructed  attractors  as  a  means  of  creating  dynamically-based  geometric  objects  that  have  a 
more  subtle  structure,  which  may  then  be  used  to  generate  state-spacc-based  predictive  models. 

Therefore,  the  idea  is  that  from  these  reconstructed  state-space  attractors,  a  baseline  predictive 
model  (generated  from  a  baseline  attractor,  X)  can  be  employed  to  make  forecasts  that  can  be 
compared  to  the  current  observed  state  (e.g.  a  test  attractor,  Y)  to  determine  if  a  change  has  taken 
place.  When  a  structure  has  been  damaged,  the  geometry  of  the  attraetor  is  presumed  to  change 
from  its  reference  state.  Thus,  one  can  use  statistically-averaged  features  characterizing  the 
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temporal  evolution  of  points  on  a  baseline  attractor  as  a  predictor  for  how  points  will  evolve  on  a 
test  attractor.  Then,  the  differences  between  the  predicted  evolution  and  the  observed  (measured) 
evolution  of  the  test  attractor  may  be  compared  in  various  ways  to  form  a  prediction  error 
feature.  Changes  to  the  prediction  error  thus  reflect  a  loss  of  dynamical  similarity  from  a 
previous  state  that  could  indicate  the  appearance  of  damage  as  it  manifests  itself  in  the 
observable  dynamic  response. 

First,  a  set  of  F  randomly  and  uniformly  selected  initial  conditions,  x„,  is  accumulated  on  a 
baseline  attractor,  X,  such  that  they  statistically  cover  the  entire  attractor  in  the  extent  of  its  state- 
space  geometry.  Next,  initial  conditions,  y„,  on  the  test  attractor,  Y,  are  determined  such  that  they 
correspond  to  the  points,  x„,  on  X.  These  points,  y„,  can  be  related  to  the  initial  conditions  on  X 
cither  spatially  (occupying  the  same  location  in  state  space)  or  temporally  (equivalent  time 
indices).  Next,  the  k  nearest  neighbors  (x^n  and  y9n .)  to  these  initial  conditions  are  found, 
separated  in  time  from  the  x„  and  y„  by  a  Theiler  window  in  order  to  minimize  artificial  near- 
time  correlations  (Theiler  1986).  Consequently,  the  neighbors  arc  purely  geometrically  correlated 
and  insensitive  to  sampling  rate.  The  size  of  the  set  of  initial  conditions  (also  called  fiducial 
points),  F,  should  be  related  to  the  total  number  of  points  on  the  attractor,  such  that  the  statistics 
generated  by  the  feature  are  insensitive  to  the  addition  of  successive  fiducial  points.  Further,  the 
number  of  k  will  depend  on  the  dynamics  and  inherent  noise  of  the  system.  Classic  choices  for  F 
and  K  arc  M10  and  Ml 000  respectively,  where  N  is  the  number  of  points  on  the  attractor 
(Pccora  and  Carroll  1996).  The  Theiler  window,  h,  is  typically  a  function  of  the  delay  choice, 
c.g.  h  =  2t,  with  the  only  stipulation  that  it  be  larger  than  T. 


Then,  the  neighbor  sets  are  time-incremented  L  time  steps  in  the  future.  The  local  prediction 
error  feature,  En(X, Y)  is  the  Euclidean  distance  between  the  evolved  fiducial  points  or 
neighborhoods  from  the  baseline  and  comparison  attractors  for  a  given  method  of  fiducial  point 
correlation  (temporal  or  spatial) 


E„(X,Y)=  kY+t-*.XJ 


(37) 


where  tj>fL  is  cither  the  initial  condition  or  centroid  of  the  neighborhood  on  the  comparison 

attractor,  evolved  L  time  steps  from  the  initial  time  n.  Likewise,  0,^  is  the  time-evolved  initial 
condition  or  centroid  of  the  neighborhood  of  the  baseline  attractor.  A  graphical  illustration  of  this 
formulation  of  this  feature  is  displayed  in  Figure  18.  The  parameter,  L,  often  referred  to  as  the 
prediction  horizon,  is  often  chosen  somewhat  heuristically  to  be  equal  to  one  or  some  value  less 
than  T.  The  feature  used  for  SHM  purposes  will  typically  be  some  global  statistic  of  the  F 
calculated  prediction  error  values,  e.g.  the  mean  and/or  standard  deviation.  This  attraetor-based 
approach  was  tried  in  the  ultrasonic  regime  in  an  application  to  detect  corrosion  in  an  aluminum 
plate,  discussed  in  the  next  subsection. 
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Baseline  Attractor,  X  _ Comparison  Attractor,  1 
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Figure  18  State-state  based  nonlinear  prediction  error. 


4.4.  Application  to  Corrosion  Detection  in  an  Aluminum  Plate 

Testing  was  performed  on  a  4'  x  4'  x  0.0625”  plate  made  of  aluminum  6061-T6.  During  testing 
the  plate  is  suspended  above  the  ground  in  a  free-free  orientation  using  wire  and  surgical  tubing. 
Nine  piczoccramic  patches  of  radius  0.125”  were  attached  to  the  plate  using  cyanoacrylate  with  a 
grid  spacing  of  12  inches.  Wires  from  each  piezoceramic  patch  were  routed  into  a  National 
Instruments  (NI)  PXI-2527  multiplexer  set  up  in  dual  16x1  2-wire  mode.  This  is  done  so  that 
both  the  actuator  and  sensor  patch  can  be  quickly  specified  through  switching  software  instead  of 
hardwiring  the  particular  pair  for  each  test.  The  excitation  signal  specified  above  is  sent  from  a 
NI  PXI-5412  signal  generator  at  a  rate  of  25  MS/s  with  14-bit  resolution  with  a  voltage  range  of 
±9  V.  This  signal  is  sent  through  a  Krohn-Hite  7602  wideband  power  amplifier  and  exits  with 
the  above  pictured  voltage  range  of  approximately  ±28  V.  The  multiplexer  routs  this  amplified 
signal  to  the  appropriate  piezoceramic  patch  and  the  ensuing  response  is  received  at  a  given 
sensing  piczoccramic  patch.  This  sensed  signal  then  goes  through  the  multiplexer  to  a  NI  PXI- 
5122  high-speed  digitizer  that  also  samples  at  25  MS/s  with  14-bit  resolution.  This  process  is 
repeated  for  all  possible  sensor  paths  with  the  multiplexer  switching,  synchronization,  and  data 
acquisition  being  controlled  by  a  custom  created  LabVIEW  virtual  instrument  file.  This 
experimental  setup  can  be  seen  in  Figure  19  below. 

Each  sensor  signal  is  averaged  10  times  and  only  the  first  1500  points  (60  microseconds)  of  the 
recorded  wave  packet  is  kept.  Only  1500  points  are  kept  so  that,  to  the  degree  it  is  possible,  only 
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the  direct  path  between  the  actuator/sensor  pair  is  interrogated.  If  more  points  were  kept  the 
signal  would  contain  more  reflections  from  indirect  paths.  This  process  is  then  repeated  100 
times  with  successive  records  added  to  the  end  of  the  last  one  so  that  a  total  record  length  of 
150,000  points  is  achieved.  Each  record  is  broken  up  into  4  runs  of  length  37,500  points.  The 
signal  is  extended  and  split  up  into  runs  in  order  to  form  more  robust  statistical  features.  All  runs 
are  normalized  by  their  standard  deviation  before  embedding  in  order  to  remove  any  possible 
environmental  variation  from  the  signal.  The  number  of  fiducial  points  used  for  each  prediction 
error  algorithm  is  10%  of  the  total  number  of  points,  or  3,750.  The  number  of  nearest  neighbors 
used  to  calculate  the  mass  centroid  is  0.1%  of  the  total  number  of  points,  or  37. 


Figure  19  Corrosion  experiment  testbed. 


Table  5:  Actuator/scnsor  pairs  for  each  interrogated  path  length 


|  Path  Length  12" 

Path  Length  17" 

Path  Length  27" 

(1— >2) 

(1 — >5) 

(1 — >6) 

(2— >3) 

(2— >4) 

(1— >8) 

(1 — >4) 

(2— >6) 

(2 — >7) 

(2 — >5) 

(3 — >5) 

(2 — >9) 

(3— >6) 

(4 — >8) 

(3— >4) 

(4— >5) 

(5— >7) 

(3 — >8) 

/— N 

V 

ON 

(5— >9) 

(4 — >9) 

(4— >7) 

(6— >8) 

(6— >7) 

oc 

A 

(6 — >9) 

(7 — >8) 

(8->9)  _  _ 

_ ■ _ i 

Table  5  is  a  summary  of  the  actuator/scnsor  pairs  that  arc  employed  in  this  experiment.  The  three 
path  lengths  used  along  with  a  corresponding  actuator/sensor  pair  are  12  inches  (1 — >2),  17 
inches  (1 — >5),  and  27  inches  (1 — >6).  Pairs  such  as  (1 — >3)  are  not  used  because  that  path  is 
effectively  covered  by  (1 — >2)  and  (2 — >3).  It  should  also  be  noted  that  for  every  listed  pair  the 
reverse  of  that  pair  is  also  used  for  redundancy,  e.g.  (1 — >2)  and  (2 — >1). 
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Three  runs  of  data  are  taken  with  the  plate  in  the  undamaged  condition  at  different  times  to 
attempt  to  account  for  environmental  variation.  The  first  run  is  used  as  the  comparison  between 
all  other  cases  and  the  second  two  runs  become  baseline  1  and  baseline  2  in  subsequent  graphs. 
Damage  is  applied  to  the  plate  through  the  process  of  electrolytic  corrosion.  In  this  experiment, 
the  aluminum  plate  acts  as  the  anode  and  a  coil  of  copper  is  used  as  the  cathode.  A  DC  voltage  of 
5V  and  1A  is  attached  to  both  the  cathode  and  the  anode  which  are  submerged  in  a  saturated 
solution  of  salt  water.  The  locations  of  the  applied  corrosion  damage  arc  shown  in  Figure  19. 
Damage  eases  1  and  2  are  applied  in  the  same  location  at  the  midpoint  between  piezoeeramie 
patches  2  and  5.  Damage  1  is  created  by  allowing  the  corrosion  process  to  take  place  for  one 
hour.  Electrolytic  corrosion  is  also  allowed  to  act  for  1  hour  in  damage  case  2.  Data  is  taken  in 
between  the  successive  application  of  corrosion  damage  at  this  location.  The  two  stages  of 
damage  can  be  seen  in  Figure  20. 


Damage  case  3  consists  of  two  separate  damage  locations.  Damage  3a  (same  size  as  damage  2)  is 
located  near  aetuator/sensor  7  and  is  created  by  corroding  the  plate  for  1  hour.  Damage  3b  (same 
size  as  damage  2)  is  located  between  patches  8  and  9  but  is  only  allowed  to  corrode  for  20 
minutes.  Thus  damage  3b  is  of  much  smaller  depth  than  all  other  damage  eases.  These  two 
damage  eases  are  pictured  in  Figure  21. 


Figure  21  Damage  case  3a  (left)  and  3b  (right). 
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Results:  12-inch  path  length 

Figure  22  depicts  the  nonlinear  state-space  prediction  error  for  all  damage  levels  as  well  as  their 
associated  10a  and  90'“  quantile  error  bars.  Redundant  paths,  such  as  (2 — >1)  for  (1 — >2),  have 
been  omitted  as  their  results  essentially  mirror  each  other  and  leaving  them  out  results  in  a 
clearer  graph.  The  solid  black  line  represents  the  path  (2 — >5)  and  shows  a  elear  monotonic 
increase  in  prediction  error  from  baseline  to  damage  1  to  damage  2.  This  is  an  excellent  result  as 
it  shows  that  this  damage  detection  algorithm  cannot  only  identify  the  existence  of  damage,  but 
that  the  nonlinear  prediction  error  will  increase  for  increasing  damage.  This  is  an  important  point 
because  if  prediction  error  did  not  have  a  positive  correlation  w'ith  damage  size  this  method 
would  have  no  potential  to  identify  the  extent  of  damage  in  any  particular  ease.  Damage  3a  is 
shown  to  be  close  enough  to  PZT  7  that  both  paths  involving  that  patch  show  a  large  increase  in 
prediction  error  for  damage  case  3.  The  preceding  damage  locations  all  have  fairly  significant 
amount  of  damage  but  damage  3b  is  a  corrosion  of  only  the  surface  of  the  plate  and  does  not 
have  much  depth.  The  path  (8 — >9)  shows  separation  from  all  other  undamaged  paths  for 
damage  ease  3  which  shows  promise  that  this  method  can  be  used  to  identify  and  loeate  very 
small  corrosion  damage.  However,  looking  at  the  two  previous  damage  eases  we  ean  see  that  the 
prediction  error  level  for  damage  3b  is  actually  less  than  level  for  several  undamaged  paths.  This 
increase  in  prediction  error  for  the  undamaged  paths  for  damage  eases  1  and  2  are  likely  not 
caused  by  the  damage  itself  because  the  prediction  error  returns  to  the  baseline  level  for  damage 
ease  3.  In  fact,  by  examining  the  time  histories  themselves  we  ean  see  that  these  differences 
appear  to  arise  due  to  extremely  slight  differences  in  the  synchronization  and  digitization 
performed  by  the  National  Instruments  PXI  data  acquisition  system.  With  this  knowledge  we 
may  be  able  to  eorreet  for  this  problem  in  subsequent  testing  but  it  remains  to  be  seen  whether 
this  increase  in  prediction  error  is  an  unavoidable  byproduct  of  the  current  test  configuration. 
Therefore,  while  damage  case  3b  has  a  prediction  error  that  falls  inside  the  range  of  several 
undamaged  paths  it  is  still  promising  that  it  is  able  to  achieve  separation  from  all  undamaged 
paths  for  the  final  measurement. 


Figure  22  12-inch  path  length  prediction  error  mean  and  90%  quantile  bars. 
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Results:  1 7-inch  path  length 

Figure  23  shows  the  prediction  error  vs.  damage  case  plot  for  the  1 7"  path  length.  Again,  damage 
3a  is  clearly  evident  in  the  path  (5 — >7).  As  one  would  expect  we  do  not  sec  any  results  in  this 
plot  for  damage  3b  because  it  is  small  and  does  not  lay  one  of  the  interrogated  paths.  For  this 
path  length  damage  case  1  yields  some  interesting  results.  We  can  see  from  Figure  20  that 
damage  1  is  nearly  symmetric  about  the  line  from  PZT  2  to  PZT  5.  However,  the  prediction  error 
results  show  damage  for  those  paths  to  the  left  of  the  damage,  (1 — >5)  and  (2 — >4),  but  no 
apparent  damage  for  the  corresponding  paths  to  the  right  of  the  damage,  (2 — >6)  and  (3 — >5). 
This  result  appears  to  be  counterintuitive  until  we  examine  a  larger  picture  of  damage  case  1  as 
seen  in  Figure  7.  This  picture  shows  several  large  dark  spots  on  the  plate  to  the  left  of  the 
corrosion  damage  that  were  caused  because  the  author  was  not  careful  to  wipe  up  excess  salt 
water  that  had  spilled  onto  the  plate  before  the  electrolytic  corrosion  process  was  started. 
Because  of  the  symmetric  nature  of  damage  1  we  must  deduce  that  the  elevation  of  prediction 
error  level  for  the  paths  (1 — >5)  and  (2 — >4)  is  a  result  of  the  slight  corrosion  that  occurred 
outside  the  intended  damage  area  due  to  the  spilled  salt  water.  This  assertion  is  backed  up  by 
examining  the  prediction  error  levels  for  damage  case  2.  If  the  aforementioned  paths  had  been 
"sensing"  the  actual  damage  1  we  would  assume  from  our  results  for  the  12"  path  length  that  a 
larger  damage  would  result  in  larger  prediction  error  for  this  case.  For  damage  case  2,  however, 
there  is  no  change  in  prediction  error  from  the  previous  damage  case  which  means  that  the  higher 
prediction  errors  for  paths  (1 — >5)  and  (2 — >4)  could  only  have  been  produced  by  the  saltwater 
corrosion  spots.  This  is  an  important  point  because  it  means  that  this  damage  detection  method 
should  be  able  to  locate  damage  better  than  it  appears  in  Figure  23.  If  the  saltwater  corrosion 
spots  were  not  present  then  damage  cases  1  and  2  would  not  be  detectable  using  the  17"  path 
length. 
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Figure  24  Saltwater  corrosion  spots  around  damage  1. 

Results:  2  7-inch  path  length 

In  Figure  25  below  we  again  see  that  damage  3a  is  easily  identified 
by  the  actuator/sensor  pair  that  includes  PZT  7  and  that  damage  3b 
is  unable  to  be  identified.  Damage  ease  1  shows  no  significant 
increase  for  either  of  the  two  paths,  (1 — >6)  and  (3 — >4),  that 
interseet  at  the  damage  location.  This  means  that  damage  1  is  too 
small  when  the  propagating  wave  has  to  travel  27"  between  the 
sensors  and  13.5"  in  either  direction  before  the  wave  interacts  with 
the  damage.  This  shows  a  limitation  for  such  long  path  lengths  but 
as  can  be  seen  for  damage  case  2,  if  the  damage  is  large  enough  a 
long  path  length  can  still  identify  the  existence  of  damage.  For  damage  case  3  we  again  sec  an 
increase  in  the  undamaged  prediction  error  level  that  partially  obseures  the  identification  of 
damage  2. 


Figure  25  27-inch  path  length  prediction  error  mean  and  90%  quantile  bars. 


Comparison  to  standard  metrics 

For  comparison,  a  five  peak  100  kHz  toneburst  excitation  was  sent  from  PZT  2  to  PZT  5  and 
PZT  8  to  PZT  9  ten  times  while  the  plate  was  in  the  undamaged  condition.  For  damage  case  1 
and  damage  case  3b  three  more  exeitation  responses  were  recorded.  Figure  26  shows  the 
averages  of  the  normalized  baseline  time  history  (solid  line)  and  the  averages  of  the  normalized 
damage  case  time  history  (dotted  line)  for  the  two  small  damage  eases.  It  is  clear  that  time-of- 
arrival  and  wave  attenuation  methods  would  be  able  to  identify  damage  easily  for  damage  ease  1 . 
However,  damage  case  3b  shows  almost  no  difference  between  the  baseline  and  damaged 
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waveforms.  This  result  shows  that  the  high  frequency  chaotic  excitation  used  with  spatio- 
temporal  prediction  error  feature  may  have  a  greater  sensitivity  to  guided  wave  metrics  because, 
as  we  noted  previously,  our  method  is  able  to  achieve  separation  between  all  undamaged  paths 
and  those  paths  that  contain  damage  ease  3b. 
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Figure  26  AO  wave  for  100  kHz  toneburst  between  sensors  2  and  5  for  damage  1 
interrogation  (left)  and  between  sensors  8  and  9  for  damage  3b  interrogation  (right). 

The  dark  line  is  the  baseline  and  the  dotted  line  is  the  damaged  case. 

5.  SUMMARY 

This  project  further  developed  and  evaluated  a  number  of  novel  tools  for  implementing  in-situ 
ultrasonic  interrogation  for  damage  detection  and  classification  in  metallic  structures  that  would 
be  of  use  to  the  US  Navy.  These  tools  included  (i)  physies-based  strueture/defcet  modeling 
using  Scmi-Analytieal  Finite  Element  (SAFE)  and  Loeal/Global  approaches  to  investigate 
specific  fcature/wave  properties;  (ii)  an  impact  detection  algorithm;  and  (iii)  data-based 
modeling  approaches  (using  time  scries  and  generalized  state-space  approaches)  for  defect 
detection,  localization,  and  assessment.  Experiments  were  conducted  in  all  eases  on  aluminum 
test  structures  to  in  notch  (simulated  crack)  detection,  bolted  joint  loss,  and  corrosion  detection. 

These  complementary  physics-based  and  data-based  tools  have  shown  great  promise  as  a 
candidate  approach  for  further  development  that  could  lead  to  prototype  systems.  Outstanding 
research  issues  that  must  be  investigated  include  (i)  robustness  to  load  and  environmental  state; 
(ii)  sensor/aetuator  bonding  normalization;  and  (iii)  sensor/actuator  architecture  optimal 
placement  strategies.  These  points  can  form  the  basis  of  a  more  thorough  proposal  that  the 
UCSD  team  will  submit  to  ONR,  given  that  this  project  was  seed-level. 
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